Stochastic cellular automata model of neural networks 
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We propose a stochastic dynamical model of noisy neural networks with complex architectures 
and discuss activation of neural networks by a stimulus, pacemakers and spontaneous activity. This 
model has a complex phase diagram with self-organized active neural states, hybrid phase transitions, 
and a rich array of behavior. We show that if spontaneous activity (noise) reaches a threshold level 
then global neural oscillations emerge. Stochastic resonance is a precursor of this dynamical phase 
transition. These oscillations are an intrinsic property of even small groups of 50 neurons. 
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I. INTRODUCTION 

Understanding the dynamics and structure of neuronal 
networks is a challenge for biologists, mathematicians 
and physicists. Neurons form complex networks of con- 
nections, where dendrites and axons extend, ramify, and 
form synaptic links between neurons. Due to long axons 
the structure of a typical neuronal network has small- 
world properties In particular, neuronal networks 
in mammalian brains have short path lengths, high clus- 
tering coefRcients, degree correlations and skewed degree 
distributions \S\ . Complex architectures of this kind are 
known to strongly influence processes taking place on 
networks [i-[7]. Complex wiring of neurons may be im- 
portant for the emergence of oscillations and synchrony in 
the brain Q. Apart from this highly heterogeneous and 
compact structure, neural networks are noisy Q. This 
makes a stochastic approach to neuronal activities un- 
avoidable d, Q . Intuitively, noise is damaging, however 
in neural networks noise can play a positive role, support- 
ing oscillations and synchrony ^, i^] or causing stochastic 
resonance [HI ■ According to experimental data, os- 
cillations and stochastic resonance may be considered as 
"noise benefits" [ll|. The origin of these phenomena, 
mechanisms and functions of oscillations in neural net- 
works are topical problems of great importance for the 
understanding of brain function [ll| . Cultured neural 
networks provide well-controlled systems for in vitro in- 
vestigations [13]. Despite their simplicity, these cultured 
networks demonstrate an extremely rich repertoire of ac- 
tivity due to interactions between hundreds to millions 
of neurons. However, at present there is no complete 
understanding of the dynamics of even these very simple 
neuronal networks. Recent investigations [HI reveal that 
global activation of living neural networks induced by a 
stimulus can be explained on the base of the concept of 
bootstrap percolation — a version of cellular automata — 
without going into details of neuron dynamics. 

In the present paper we propose a stochastic cellular 
automata model of noisy neural networks. Based on ex- 
perimental data we assume that activation processes are 
stochastic, i.e., neurons can be activated with a certain 
probability either by an external stimulus, spontaneously. 



or by fluctuating inputs from active presynaptic neurons. 
These networks include two neural populations, excita- 
tory and inhibitory neurons, and have a complex network 
architecture, i.e., the small world property and hetero- 
geneity are taken into account. We consider model neu- 
rons which fire regular trains of spikes with a constant 
frequency. The stochastic dynamics of these networks 
takes into account processes of spontaneous neural ac- 
tivity, which plays the role of noise, the activation of 
neurons by a stimulus or neural pacemakers, and inter- 
actions between neurons. With this model we aim to 
understand the role of noise in the emergence of oscilla- 
tions and the origin of stochastic resonance. Although 
the model is simple, it demonstrates various patterns of 
self-organization of neural networks, hybrid phase tran- 
sitions, hysteresis phenomena, neural avalanches and a 
rich set of dynamical phenomena driven by noise: decay- 
ing and stable oscillations, and stochastic resonance. 

Using exact analytical methods and simulations of the 
stochastic dynamics of this model, we demonstrate that 
noise can play a constructive role in neural networks. We 
show that at a critical level of noise a neural network un- 
dergoes a dynamical phase transition from a state with 
incoherent neurons to a state with synchronized neurons 
and global oscillations. Oscillations of neural populations 
emerge if spontaneous neural activity (noise) is above a 
critical level. Stochastic resonance is a precursor of global 
oscillations. At a given spontaneous neural activity, a 
critical fraction of neural pacemakers can also stimulate 
oscillations. We consider several mechanisms leading to 
global oscillations in neural populations: the difference 
in dynamics of excitatory and inhibitory neurons or the 
existence of synaptic delays. These mechanisms lead to 
similar oscillations. We also show that global oscillations 
are intrinsic properties of the neural networks under con- 
sideration. One should note that these oscillations are 
nonlinear waves with a certain amplitude and a specific 
shape which are determined by the structural and dy- 
namical parameters. They do not depends on initial con- 
ditions in contrast to waves in linear dynamical systems. 
We demonstrate that the network structure plays an im- 
portant role. In neural networks having the structure of 
classical random networks the larger the connectivity the 
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broader is the region with global oscillations. Our sim- 
ulations reveal that oscillations are an intrinsic property 
of even small groups of neurons. 50-1000 neurons dis- 
play oscillations similar to infinitely large networks de- 
spite stochastic fluctuations which are usually strong in 
small networks. The proposed model also explains a dis- 
continuous transition in the activation processes of living 
neural networks observed experimentally in |12| . Neural 



avalanches precede this transition. Simulations support 
our analytical solution. 



II. MODEL 



a) 



b) 



FIG. 1: Firing rate v versus input V: (a) type 1 neuron; (b) 
type 2 neuron; (c) the step function approximation used in 
the present paper. 



Neurons demonstrate various types of spiking behav- 
ior in response to a stimulus at firing threshold, see, for 
example, [13h15| . Type 1 neurons show a continuous 
transition from an inactive state to an active state with 
an arbitrary low firing rate when the input current is 
above a threshold input (see Fig. [1^). For example, cor- 
tical excitatory pyramidal neurons exhibit this behavior. 
Frequencies of tonic spiking of type 1 neurons lie in the 
range from 2 Hz to 200 Hz, or can be even higher than 
200 Hz. The maximum firing rate is set by the refractory 
period of a neuron. Type 2 neurons show a discontinu- 
ous transition to a nonzero firing rate above a threshold 
input (see Fig. [Dd). They fire in a relatively narrow 
frequency band. For example, Hodgkin-Huxley neurons 
demonstrate type 2 neural excitability. Type 2 neurons 
fire spikes with frequency about 40 Hz and higher. Fast- 
spiking inhibitory interneurons in the rat somatosensory 
cortex fire in the frequency range 20-61 Hz p^ . Neurons 
with type 2 dynamical behavior may play an important 
role in synchronization of neural activity [13] • Several 
models have been proposed to describe the dynamics of 
individual neurons (see, for example, [Tsl Il8l - l20| ). 

In the present paper we only consider regular spiking 
neurons. We approximate the frequency-current response 
by the step function (see Fig. [TJ;). Active excitatory 
and inhibitory neurons fire trains of spikes with a con- 
stant frequency v which is the same for all neurons and 
does not depend on the input. \i tv > \ then during 
an integration time r (the membrane time constant) a 
postsynaptic neuron receives \rv\ spikes from an active 
presynaptic neuron, where \A\ stands for the integer part 
of a number A. It is assumed that the spike duration 
(about I ms) is much smaller than r. The membrane 
time constant r can range from 1 to 100 ms [2l|. For ex- 
ample, for a typical integration time r = 10 ms we must 
have V >100 Hz. 

Let us consider a neural network with two types of 
neurons: excitatory and inhibitory neurons (see below). 
The total number of neurons is N . The fractions of ex- 
citatory and inhibitory neurons are and gi = 1 — ge, 
respectively. Neurons are linked by directed edges and 
form a network with an adjacency matrix anm where 
n,m = l,2,...,iV. An entry a„m is equal to 1 if there 
is an edge directed from neuron n to neuron m, other- 



wise a„„i — 0. Each neuron can be in either an active 
or inactive state. Active neurons fire regular trains of 
spikes, as discussed above. We assume that there is no 
phase correlation between trains of spikes generated by 
different neurons. We define s„(<) = 1 if neuron n is ac- 
tive at moment t, and Sn{t) = if this neuron is inactive. 
In our model, these binary variables play an auxiliary 
role. During the integration time r a postsynaptic neu- 
ron receives and integrates spikes from active presynaptic 
excitatory and inhibitory neurons. First we consider the 
case TV > 1. The total input Vn{t) (post-synaptic po- 
tential) at neuron n is the sum of inputs from nearest 
neighbor (presynaptic) neurons: 

^vnnJrani (1) 

where synaptic efficacy J^n = ±J if neuron m is exci- 
tatory or inhibitory, respectively. We assume that all 
synapses of excitatory neurons are excitatory, and all 
synapses of inhibitory neurons are inhibitory. This is the 
so called Dale's principle [13]. Recently, the importance 
of Dale's principle for dynamics and pairwise correlations 
in neural networks was discussed by Kriener et al (23j . In 
our model, the dynamics do not change qualitatively if 
\Jmn\ are different for these two populations of neurons. 
Note however that there are physiological reasons for the 
fact that the magnitudes of inhibitory efficacies are usu- 
ally larger than excitatory efficacies (see, for example, 
[2l|). Active excitatory (inhibitory) presynaptic neurons 
give positive (negative) inputs to a postsynaptic neuron, 
while inactive neurons give no input. For example, the 
input from k active excitatory and I inhibitory neurons 
\sV — [Tv\Jk — [tv\JI. We suppose that this input acti- 
vates the postsynaptic neuron if V is at least a threshold 
value Vth- This gives the following condition: 

k -l>VL = Vthl[Tv]J (2) 

which we will use below. Notice that £7 is a dimension- 
less parameter. The dimensionless threshold 51 is of the 
order of 15-30 in living neural networks and about 
30 — 400 in the brain. Even if biological neurons have a 
variable threshold 14], for simplicity we assume that the 
threshold does not depend on the prior activity. 
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In our stochastic model we assume that the states of 
neurons at each moment t are determined by the follow- 
ing rules: 

(i) An excitatory (inhibitory) neuron is activated at a 
rate /e {fi) either by a stimulus or spontaneously 
(spontaneous activity). 

(ii) In addition, an excitatory (inhibitory) neuron is ac- 
tivated at a rate /zie (fJ-u) by nearest neighbor ac- 
tive neurons if the total input V{t) at this neuron 
is at least a threshold value Vth, i.e., V{t) > Vth- 

(iii) An activated excitatory (inhibitory) neuron is in- 
activated (i.e., it stops firing) at a rate /iie {^J■li) if 
the total input V{t) becomes smaller than Vth- 

(iv) An activated excitatory (inhibitory) neuron spon- 
taneously stops firing at rate fi2e {p-2i)- 

In the brain, neurons receive fiuctuating inputs and 
generate spike trains [8]. We represent the activation 
by fluctuating inputs as the stochastic process (ii) with 
the rates and fin which can be of the order of the 
average firing rate. This determines the time scale in 
the model. Even if the total input is on average larger 
than Vth, it sometimes falls below Vth- As a result, the 
neuron stops firing. Process (iv) is meant to represent 
this. The biophysical meaning of the model parameters, 
assumptions and approximations which are the basis of 
our model, are discussed in Sec. IVIII For other models 
with binary variables see (23-[27i and in the review [20| . 

In order to describe the dynamics of neural networks, 
we introduce a probability pii^ (t) that neuron n of type a 

is active at time t. Let us define the mean values of (t) 
for excitatory, a — e, and inhibitory, a = i, populations: 

Pait)^Y.Pn\t)/i9aNl (3) 
n 

where the sum is over neurons of type a, ga is their frac- 
tion. We name Pe{t) and pi{t) "activities" of the exci- 
tatory and inhibitory populations. On the other hand, 
Peit) and pi(t) are the respective probabilities that a ran- 
domly chosen excitatory or inhibitory neuron is active at 
time t. We consider neural networks whose structure is 
of a sparse random uncorrelated directed network. These 
networks are small worlds and can have an arbitrary de- 
gree distribution. They are often considered as a good 
approximation to real networks 2]. The advantage of 
these networks is that they can be studied analytically 
by use of mean-field theory and easily modeled for sim- 
ulations. However, they do not take into account the 
high clustering coefficient and degree correlations of real 
neural networks Q. Though the mean- field approach is 
based on the tree-like approximation, it takes into ac- 
count exactly the heterogeneity of networks and large 
feedback loops Q. 



III. BASIC RATE EQUATIONS 

Let us derive dynamical equations for the activi- 
ties pe{t) and Pi{t). We introduce the probabilities 

eiPeit), Pi{t)) and i{pf.{t) , pi{t)) that at time t the to- 
tal input to a randomly chosen excitatory or inhibitory 
neuron, respectively, is at least $7. If at time t an excita- 
tory neuron is inactive, which takes place with probabil- 
ity 1 — Pe{t), then an external field activates this neuron 
at a rate /g. This gives a contribution 

fe[l-Pe{t)] (4) 

to the rate Pe{t) = dpe{t)/dt. If at time t the total input 
to an inactive neuron is at least il, which takes place with 
probability "i! e{peit), pi(t)), then this neuron is activated 
at the rate /iie. This gives one more positive contribution 

flle[l- Pe{me{Pe{t),p.M)- (5) 

If at time t an excitatory neuron is active, which 
takes place with probability Pe{t), and the total in- 
put from activated nearest neighbor excitatory neurons 
is smaller than £7, which takes place with probability 
1 — ^'e(/Ce(i), /3i(i))j then such an active neuron becomes 
inactive at the rate pie- The active neurons also can stop 
spontaneously firing with rate p2e- These processes give 
two negative contributions: 

- piePe{t)[l - *e(Pe(t), P^{t))] - fl2ePe(t)- (6) 

Summing all contributions, we obtain a rate equation, 

Pa{t) = fa - VaPa{t) + P'la'^ a{Pe{t) , Pi{t)) . (7) 

Here Va = fa+P-ia+P-2a, and a = e,i. 

To clarify the relative role of activation and deactiva- 
tion processes, we rewrite Eq. ([7]) as follows: 

Pa/Va=Fa{l-Qa)-Pa + {^-Fa){l-Qa)^a{Pe, P^), (8) 

where pa — pa{t)- The dimensionless parameters Fa = 
falifa + Mia) and Qa = P2a/va determine the relative 
strength of stimulation and the spontaneous deactivation 
of neurons. The rates Ve and Vi set the time scale. 

The probabilities 'I'e and VP^ are determined by the net- 
work structure. Below we will study a directed classical 
random graph which is the simplest and representative 
model of sparse uncorrelated complex networks [2, 
These random graphs share the properties of sparse un- 
correlated random networks with a finite second moment 
of the degree distribution. They are small worlds and 
have a mean shortest distance which increases as the log- 
arithm of the number of vertices, in contrast to a three 
dimensional system where a mean shortest distance in- 
creases as the cube root of the size. Due to simplicity, 
classical random graphs are often used to study dynamics 
of systems having a complex network structure ■ 
In contrast to real networks, sparse random uncorrelated 
networks and in particular classical random graphs have 
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zero clustering coefScient due to their tree-like structure 
and negligible (in some cases, weak) degree-degree cor- 
relations between neighboring nodes in the infinite size 
limit. Understanding the strength of the clustering and 
degree correlations on dynamics of systems with complex 
network architecture is an open problem in the theory of 
complex networks Recent investigations of vari- 

ous dynamical models on complex networks show that 
in many cases networks with clustering demonstrate dy- 
namics qualitatively similar to tree-like networks. In 
many cases degree-degree correlations also do not quali- 
tatively change the dynamics. This challenging problem 
is discussed in detail in the recent review [5j. 

In a classical random graph a directed edge between 
each pair of N neurons is present with a given probability 
c/N. The parameter c is the mean input and output 
degrees. The probability i?„(c) that a neuron has n input 
edges is given by the binomial distribution: 

i?n(c)-C^(^r(l-^r-" (9) 

where — Nl/{N — n)\nl is the binomial coefficient. 
We will study analytically large networks with N ^ 1. In 
the infinite size limit, N — > oo, the binomial distribution 
Bn{c) approaches the Poisson distribution Pn(c), 

Pn{c) = c"e-7n!, (10) 

which is more convenient for calculations. The probabil- 
ity that a randomly chosen neuron has k active presynap- 
tic excitatory and / active presynaptic inhibitory neurons 
is Pk{gePec)Pi{giPic) ■ Heucc, in the case tv > 1, we get 

fe-O 

*e(Pe,Pi) = *j(Pe,Pi)= ^ ^ Pki9ePeC) Pl{giPiC) 
k>n 1=0 

fc>o ^ ' 

where T{k, x) is the upper incomplete gamma function 
and 51 is defined by Eq. ([2]). Notice that in the case 
of classical random graphs we have used the fact that 
there are no correlations between the number of input 
and output edges. 

In the case tv < 1, during the integration time r a 
postsynaptic neuron receives only one spike or none from 
an active presynaptic neuron. If the phase of a train of 
spikes is uncertain then all we can say is that during the 
time interval r with probability tv a postsynaptic neuron 
receives a spike from an active presynaptic neuron. In 
turn, the probability that there is no spike is 1 — tv. 
Let us assume that there is no phase correlation between 
regular spiking neurons. This is a common assumption 
at low activity rates [l^. The probability that during 
time T a neuron receives k spikes from uncorrelated n 
regular spiking presynaptic neurons is 

C^{Tv)\l-Tvr-\ (12) 



In the case of a classical random graph, the probabil- 
ity that during the integration time r a randomly cho- 
sen neuron receives k spikes from active excitatory or 
inhibitory neurons is given by the Poisson distribution: 

oo 

J2Pn{9aPaC)Cj:{Tv)'^{l~TVr-'^ = P.igaPaTVc), (13) 
n— /e 

where a = e,i for excitatory and inhibitory neurons, re- 
spectively, k spikes from excitatory and / spikes from 
inhibitory neurons activate a postsynaptic neuron iiV — 
Jk — Jl > Vth- Using the probability Eq. ([T^. one can 
show that the function 5'a(pe(i), pi{t)) in Eq. ([7]) is given 
by Eq. (1111) if the mean degree c is replaced by tvc, and 
a threshold ft — Vth/ J is used. Therefore, the effective 
mean input degree is decreased while the effective thresh- 
old is increased in comparison to the case tv > 1. Note 
that if trains of spikes generated by presynaptic neurons 
are correlated, then Eq. (IT^ is invalid. Spikes acting in 
concert can activate a postsynaptic neuron more effec- 
tively. 

One can use another approach. The stochastic rules 
(i)-(iv) lead to a rate equation for the activity p'^\t) of 
single neuron n of type a with g„ = Omn presynaptic 
neurons for a given adjacency matrix anm'- 

p^:\t)^fa-vap\:\t) 

t)], (14) 

{s„=0,l} rii 

where Vn = [Tv]^^SmO-mnJmn is the input at neuron 
n from presynaptic neurons m at [tv] > 1, Q{x) is the 
Heaviside step function, pm{s„i—0,t) = 1 — p[n {t) and 

Pm{sm=^,t) = Pm\t) are the probabilities that presy- 
naptic neuron m is inactive or active at time t, respec- 
tively. The last term in Eq. ([T^ is the probability that 
the input at neuron n is at least the local threshold Vth {n) 
at time t. These equations describe a neural network 
with a given adjacency matrix anm, arbitrary synaptic 
efhcacies Jnm and arbitrary local thresholds Vthiji). In 
the case of the classical random graph in the infinite size 
limit, for the uniform case \Jnm\ = 1 and Vth{n) = Vth, 
the set of TV coupled nonlinear rate equations can 
be reduced to two coupled equations for the averaged 
activities pe and pi. Summing over n in Eq. (|14p and av- 
eraging over the network ensemble, we arrive at Eqs. ([7]) 
and pT|) . We believe that the mean- field equation ([7]) 
is exact for sparse uncorrelated directed networks in the 
limit N ^ oo. Our simulations of the model on classical 
random graphs support this. Similar rate equations were 
derived for disease spreading and contact processes on 
complex networks ^2^, ^SflJ . 

Neural networks can also be activated by pacemak- 
ers (neurons that permanently fire). Let excitatory and 
inhibitory pacemakers be chosen with given probabilities 
Fe and Fi from excitatory and inhibitory neurons, respec- 
tively. The stochastic dynamics of remaining neurons 
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FIG. 2: Function 'i/a{pe,Pi), Eq. (fTT|) . versus QePe at giPi—O, 
0.3, 0.6, and 0.8. Other parameters: c = 100, D. = 20. 



(activities Pe{t) and Pi{t)) are governed by rules (ii)-(iv). 
In the same way as for Eq. ([8]), we obtain 

Pa/l^a = Fa-pa + {l-Fa){l-Qa)'i>a{Pc,Pi), (15) 

where we define pa=Fa+{l~Fa)pa{t), the total activity 
of the neural population a, a = e, i. Equations ^ and 
(|15|1 differ only by the first term on the right-hand side. 
Thus, activation by a stimulus or randomly chosen pace- 
makers produce similar effects. A similar equation at 
Qa ~ was derived using another approach in ^26j . 

In our model one can also take into account synaptic 
delays. Introduce time Tab for the transmission of a nerve 
signal from a neuron of type a to a nearest neighbor neu- 
ron of type b, where a,b = e,i. Then, in Eq. (|5]), replace 

^aiPeit),P^{t)) by a[Peit - Tea),P^{t - Tia)]. VarioUS 

sources of delays in the nervous system and their role in 
dynamics of neural networks were recently discussed by 
Ermentrout and Ko [3ll |. 

The rate equations ([7]) look similar to the rate equa- 
tions derived in the pioneer works of Wilson and Cowan 
[3^ [33! who considered the dynamics of neural popula- 
tions with excitatory and inhibitory interactions. How- 
ever, there are important differences between our model 
and the Wilson-Cowan model. Our model of interact- 
ing excitatory and inhibitory neurons is based on the 
stochastic rules of activation and inactivation of individ- 
ual neurons (these are the rules (i)-(iv) in Sec. |ll| in 
contrast to the deterministic phenomenological model in 
[H,!!!]. Using these rules, we derived the self-consistent 
rate equations ([7]). Furthermore, Wilson and Cowan 
used as relevant variables the fractions of excitatory and 
inhibitory neurons which become active per unit time. 
Within our notations these are gePe and giPi, respec- 
tively. In our approach in the case of classical random 
graphs, the fractions of active excitatory and inhibitory 
neurons, i.e., gepe and giPi, are the relevant variables. 
Also, on the base of experimental studies, Wilson and 
Cowan postulated that the subpopulation response func- 
tions have a sigmoid form. They used the standard mean 
field theory which neglects the spatial heterogeneity, and 



FIG. 3: Activity pe of excitatory neurons versus the activation 
parameter F at different fractions of inhibitory neurons gt 
from numerical solution of Eq. ((Sj at c = 20, Q = 3. The 
jump and hysteresis disappear if gi > g* ~ 0.43. Arrows 
show increasing and decreasing F. The insert shows results 
at c = 1000, Q, — 30. Our simulations confirm these results. 



assumed that all neurons are subjected to the same aver- 
age excitation of excitatory and inhibitory populations. 
In our model, the functions ^'e(pej/3i) and ^i(pe,pi) in 
Eqs. ([7]) play the role of the response functions. We cal- 
culated these functions exactly, taking into account the 
heterogeneity of the classical random graph. According 
to Eq. (jlip . these functions have a sigmoid form with 
one inflection point as a function of the parameter gePe 
in a wide range of giPi, see Fig. O One can expect a 
multimodal functional dependence with several inflection 
points if there arc several neural populations with differ- 
ent thresholds Vth- Finally, in our stochastic approach, 
the set of Eqs. (fH|) permits the study of the dynamics 
of individual neurons while Eqs. ([7]) describe the global 
activity of the neural populations. The Wilson-Cowan 
model only describes the global activity of the neural 
populations. Below we will show that the stochastic 
model as well as the Wilson-Cowan model reveal hystere- 
sis phenomena, decaying and stable oscillations in neural 
activity. 



IV. STEADY STATES AND AVALANCHES 

The steady states of the model are determined by 
Eq. ^ at Pa = 0. The steady solutions of Eq. ^ gen- 
eralize the standard bootstrap percolation to a directed 
random graph with two types of vertices. A particular 
case with gj = 0, Fg — Fi, and Qe = Qi = was studied 
in Refs. Activation processes are shown in Fig. [3] at 
F = Fe=Fi, Qe=Qi=0 when Pf, = pi. One can see that 
by increasing the activation parameter F, the activity pe 
(and Pi) undergoes a jump at a critical point Fc. A sim- 
ilar jump was observed in living neural networks in vitro 
[T^ . If F approaches Fc from below, then 

Pa = pi'^-A{Fc^Fy/\ (16) 
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where A is a coefRcient. This singular behavior evidences 
the existence of long-range correlations between neurons 
and the emergence of neural avalanches: the activation or 
deactivation of one neuron triggers the activation or de- 
activation of a large cluster of neurons. This phenomenon 
is similar to one that was found near the point of emer- 
gence of a giant fc-core [s^] . Thus the transition at Fc is a 
hybrid phase transition (one which combines a jump and 
a singularity). At F = Fc the probability G{s) that an 
avalanche has a size s, including the activating neuron, 
is 



(17) 



Similar neuronal avalanches were observed in the cortex 
[13, fsl]. Using the approach from we calculated 
G{s) exactly at 5^ = and F < F^. 



G{s) = ^!!£!fllle-"-^ 



(18) 



where the average number of inactive subcriti- 

cal postsynaptic neurons of an inactive presynaptic ex- 
citatory neuron. By definition, a subcritical neuron 
has exactly 0,-1 active presynaptic excitatory neurons. 
Successive activation of these subcritical neurons form- 
ing finite clusters leads to avalanches. We found that 
Ucr = (1 — F)d^eiPe,(^)/dpe < 1 whcrc Pe is the neural 
activity in the steady state at a given F. At the critical 
point F — Fc we have Ucr = 1- This leads to Eq. (flT)) 
which we believe is also valid for gi ^ 0. 

With increasing gi the size of the jump decreases. 
There is a special critical point g* at which the jump 
is zero, and the phase transition is continuous. There is 
no phase transition if gi > g* , or if f2 is larger than a 
critical threshold (see Fig. In Fig. [3] we display nu- 
merical results for large mean degree c=1000 and large 
il=30, and for small mean degree c=20 and small f2=3. 
Qualitatively the behavior is the same. There is a range 
of gi in which the system demonstrates bistability while 
the upper metastable state has activity pe not close to 
1 (see small hysteresis loops in Fig. [S]). However, this 
region becomes smaller in the case of large c. This indi- 
cates that with increasing c and f2, this bistability region 
decreases rapidly. Thus the hysteresis behavior crucially 
depends on having finite values of c and fl. In biological 
systems the efficacy of inhibitory synapses is larger than 
that of the excitatory ones. In our analysis we assumed 
that they are equal. Our calculations show that an in- 
crease in magnitude of the inhibitory efficacy moves the 
fraction gi of inhibitory neurons at which the interesting 
bistability region takes place into a region of biologically 
plausible values, namely about 0.2. 



V. RELAXATION AND OSCILLATIONS 

Let us consider the relaxation of neural networks to a 
steady state. We represent Pait) as pa + Spa{t) where 



5pa(t)l Pa ^ 1, and pa is the equilibrium activity of pop- 
ulation a. Linearization of Eqs. (jS]) with respect to (5 /9a(0 
gives two coupled linear equations: 



dSpa{t) 

Undt 



-5 Pa {t)+DaeSpe {t)+Dai5pi (t) , 



(19) 



where Dab=il-Fa)ii-Qa)d^ aiPe, Pi) /dpt for a, 6 = e,i. 
We look for a solution in the form Spa{t) = Aae~'^* with 
unknown Aa and 7. The solution exists if the determi- 
nant of this set of equations is zero. This condition gives 
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= i^c{Bi+B2±[{B,-B2f+4aDcM'^^}/^, (20) 



where a = Vijvc^ B\ = 1— Dee, ^2 = ct(l— i'ii)- Equa- 
tion (l20t is valid in the general case 4'e 7^ ^i- For the 
classical random graph, using Eq. ([TT|). one can prove 
that Dee, Die > while Dci,Dii < 0. Therefore 7 in 
Eq. (1201) may be a complex number in certain ranges of 
parameters c, g, F, and a. Where 1017 = 0, relaxation 
is exponentially fast with the rate 7. For example, at 
a = 1, we have 7 = vd^—Dee—Da) > 0. In this case 
7 tends to if F Fc from below as at a continuous 
phase transition. However 7 is always finite above the 
critical point Fc. If Re7 > and Im7 ^ 0, then relax- 
ation is in the form of decaying oscillations. If Re7 < 
and Im7 ^ 0, then any small deviation from a steady 
state leads to oscillations around the state with an in- 
creasing amplitude. However, in this case the linear ap- 
proximation, Eq. (1191) , is not valid, and it is necessary to 
solve Eqs. These three regions are shown in Fig. |4l 
We solved Eqs. ([8]) numerically in the case Fc = Fi — F, 
Qe = Qi = 0. We found that there is a region oi gi, which 
includes the special point g*, where Re7 < and Im7 7^ 
if 0<a<ac2=(£'ee — 1)/(1— -Dii) < 1, i.e., when inhibitory 
neurons have slower dynamics compared to the dynam- 
ics of excitatory neurons. It turns out that in this region 
the neural system displays stable oscillations around the 
steady state. Figure |4] shows that the larger the mean de- 
gree c and the threshold fi the broader is the region with 
oscillations. We obtained similar results for the model 
with synaptic delays. In particular, there is a region of 
gi where oscillations emerge at a=l and Tee=Tei=0 if 
Tic=Tii > T where T is a threshold. The firing rate pi 
in human brains is typically in the range 1 - 400 Hz. In 
our model the frequency of oscillations uJq is several times 
smaller than /ii. This gives Wo h^ the range of the waves 
observed in brain, i.e., Wq < 100 Hz. 

Replacing fa by fa{t)^fa+^aSm{ujt) in Eq. d?]), we 
study the response of the model, pa+Apa sm{ujt+(pa), to 
a small periodic stimulation, Aa fa- If F approaches 
the boundary between regions (II) and (III), see Fig. 21 
the response 



(ApaMa)Vl/[(a;-Im7)2 + (Re7)^ 



(21) 



is enhanced because Re7=0 at the boundary. There- 
fore the transition from a state with incoherent neurons 
to a state with global oscillations is a dynamical phase 
transition with a sharp boundary (in the thermodynamic 
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FIG. 4: There are three regions on the a — F plane: (I) with 
exponential relaxation; (II) with decaying oscillations; (III) 
with stable oscillations. The boundaries Qci and ac2, given 
by equations Im7 — and Re7 — 0, are shown at gi < g* 
(solid lines), and Qi > g* (dashed lines), (a) c = 20, f2 = 3, 
=0.4 and 0.47. (b) c = 1000, fl ^ 30 gi =0.475 and 0.478. 




FIG. 5: (color online). Fractions Re = Pege and Ri = pigi of 
active excitatory and inhibitory neurons versus time, (a)-(c); 
a = 1 (region (I)), (d)-(f): a = 0.4 (region (II) ). (g)-(i); 
a = 0.05 (region (III)). Solid (dashed) lines show theoretical 
Re (Ri) from Eqs. Blue (red) symbols refer to Re (Ri) 
from simulations at A*' = 10000 (1st row), 1000 (2nd row), and 
50 (third row). F = 0.05, gi = 0.4, c = 20, f7 = 3, ^2 = 0. 

limit). In our model the stochastic neural activity plays 
the role of noise while interactions between neurons pro- 
duce non-linear effects. Thus the observed strong en- 
hancement of the response is actually stochastic reso- 
nance [lll,|33]. 

VI. SIMULATIONS 

Our simulations supported the theoretical results. 
Random networks with N neurons were constructed by 



establishing directed links between any neuron i and neu- 
ron J, with probability c/N. In the initial configuration 
all neurons were inactive. The state of each neuron is 
then updated every At time units (parallel update) ac- 
cording to stochastic rules (i)-(iv). [Any other initial con- 
figuration may be also used.] The value of the time step 
At was chosen such that the probabilities fAt, fiiAt, and 
li2At of the stochastic processes (i)-(iv) in Sec. |lI]for ex- 
citatory and inhibitory neurons were sufficiently small. 
Reliable results were obtained when these probabilities 
were about 0.1 or smaller. Figure [5] represents typical 
numerical results obtained for systems of different sizes. 
All parameters used in simulations are presented in the 
caption to Fig. 5. For a given number of neurons N 
we constructed several realizations of networks, and then 
we simulated their stochastic dynamics, using the rules 
(i)-(iv) in Sec. [TTl As one would expect, for the consid- 
ered stochastic model, different runs and different realiza- 
tions of neural networks differ slightly one from another. 
With increasing N these differences become smaller and 
smaller, so these are standard run-to-run and realization- 
to-realization variations. 

Figure [S] shows a full set of regimes. One can see that 
in regimes with exponential relaxation and decaying os- 
cillations the irregular activity of neurons decreases with 
increasing N. Already at = 1000, a stimulation with 
F > Fc activates a finite fraction of neurons in agreement 
with the theory, though there are strong irregular fluc- 
tuations around the steady state. In a small network of 
50 neurons stochastic effects are strong and suppress the 
global activation. In Fig. [5] we also compare oscillations 
predicted by Eq. (|5]) to our simulations. Interestingly, 
these oscillations have a saw-tooth shape. Their period 
and shape depend on the parameters of the model such 
as F, a, c, fl, and gi. The theory and simulations are in 
very good agreement at A^ = 10000. Actually we found 
good agreement with only A^ = 1000. Surprisingly, the 
predicted oscillations emerge even in small groups of 50 
neurons where strong stochastic effects and non negligi- 
ble clustering could be expected. For c = 20 and A^ = 50 
the mean clustering coefficient is C = c/N = 0.4 0, [3, 
which is close to the value C = 0.53 found in the macaque 
visual cortex Q . This intrinsic property of small groups 
of neurons to oscillate may be very important for under- 
standing communication between neuronal groups in the 
brain [38i] . 



VII. DISCUSSION 

First let us discuss the assumptions and approxima- 
tions which are the basis of our stochastic approach to 
noisy neural networks, and explain the biological mean- 
ing of the model parameters from the point of view of 
experimental and theoretical neuroscience. 

In our model, activation of neurons by stimulus is a 
stochastic process with a characteristic time which is 
equal to the reciprocal rate In the brain, stochas- 
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ticity in activation of neurons by stimulus may appear in 
trial-to-trial variability of the first-spike latency of neu- 
rons. The first-spike latency of a given neuron is defined 
as the time from the onset of a stimulus to the time of 
appearance of the first spike. The first-spike latency can 
depend on many parameters. For example, for auditory 
neurons it depends on the amplitude and frequency of 
stimulus [22]. We suppose that the reciprocal rate l/^i 
is of the order of the mean first-spike latency of neurons. 
For simplicity, we assume that l/^i is constant and does 
not depend on the input. The first-spike latency may 
be of the order of the period of tonic spiking or much 
larger if the input is near the threshold. In the mam- 
malian cortex the latency of regular spiking neurons for 
a superthreshold input can be in the tens of milliseconds. 
This gives /^i - 10 - 400 Hz. 

In our approach it is assumed that each neuron may 
be active spontaneously, that is, it may discharge with- 
out experimenter-controlled stimulations. At the present 
time, the mechanisms and functional significance of spon- 
taneous neural activity are not well understood and it is 
a topical problem in experimental and theoretical neu- 
roscience [1, 0, [H, ESI ■ The typical spontaneous back- 
ground activity observed in the cortex is 1-5 spikes/s. 
Interactions between neurons play an important role in 
this activity (28j . Spontaneous activity in the brain may 
be mediated by intrinsic, intracellular, generated activity 
and circuit feedback mechanisms. Neural activity in one 
region of the brain may propagate to other regions, cir- 
culating in recurrent loops. For example, neurons in the 
thalamus and the cerebral cortex form recurrent loops 
(4l| . The study of spontaneous activity in neocortical 
slices [i^ gives evidence that supports both mechanisms. 
For our model a mechanism of spontaneous activity is 
unimportant. One can assume that spontaneous activ- 
ity takes place by the intrinsic mechanism. Alternatively 
one can consider the neural network as part of a large 
system from which neurons receive random inputs. In 
real neural networks only a fraction of neurons are spon- 
taneously active [ilj. In the present paper we study the 
case in which all excitatory and inhibitory neurons may 
be spontaneously active. One can show that if only some 
fraction of neurons is spontaneously active, the dynamics 
of the neural networks would be qualitatively the same. 

Furthermore, we considered the activation of neurons 
by an external stimulus as a stochastic process. Exper- 
imental work supports this assumption. For example, it 
was revealed that a moving whisker can have only a 15% 
chance of generating spikes in a neuron in the mouse 
somatosensory cortex |43i |. Unfortunately, much less is 
known about stochastic processes of spontaneous deacti- 
vation of neurons. In our model, the reciprocal rate l//i2 
is the characteristic time at which neurons stop firing due 
to irregular fluctuations on the input or due to random 
processes taking place inside the cells. Recently it was 
shown that spontaneous activity of single neurons may 
be driven by noise which can not only activate but can 
also inhibit spiking activity of neurons of both type 1 and 



type 2 gi-Iii]. 

The rates of the stochastic processes discussed above 
can be found from statistical analysis of activation and 
inactivation events in neural networks. They can also be 
measured by use of the patch-clamp technique: one can 
stimulate presynaptic excitatory and inhibitory neurons 
and then measure the probability of activation of a post- 
synaptic neuron through the distribution of first-spike 
timing times. 

The proposed stochastic model is not restricted to reg- 
ular spiking neurons. One can also analytically study 
noisy neural networks with neurons which generate ran- 
dom spike trains, for example, Poisson spike trains as 
found in recordings from neurons in vivo and in vitro [9] . 
The proposed stochastic approach can also be generalized 
to study analytically neural networks with neurons hav- 
ing the type 1 and 2 dynamical behavior shown in Fig. [1] 
for the case when correlations between presynaptic neu- 
rons may be neglected. However these generalizations 
are out of the scope of the present paper. 

We found that even a small group of neurons re- 
veals intrinsic oscillations which are robust against strong 
stochastic fiuctuations (see Fig. [SJ . It means that despite 
noise, neurons in a small group can synchronize their 
dynamics. We believe that this result opens interesting 
possibilities to study and model communication between 
different groups of neurons and the transmission of activ- 
ity from one group of neurons to another. In the recent 
review (38| . "neuronal communication between neuronal 
groups through neuronal coherence" was considered as 
a mechanism for cognitive dynamics. On the basis of 
neurophysiological data. Fries suggested that cohere ntly 
oscillating neuronal groups can interact effectively [38| . 
This idea is based on an assumption that activated neu- 
ronal groups have an intrinsic tendency to oscillate. Our 
model supports this assumption, showing that oscilla- 
tions indeed are an intrinsic property and robust against 
noise. In our model one can model communication be- 
tween neural populations or neural modules. Synchro- 
nization of neurons or groups or modules of neurons in 
the regime with oscillations can play an important role in 
this communication. Indeed, our preliminary simulation 
of interacting neural communities reveals complex pat- 
terns of neural activities. On the basis of our stochastic 
model one could study the computational role of network 
oscillations and how oscillations contribute to the repre- 
sentation of information 47] . 

Real neural networks have a scale- free degree distribu- 
tion rather than a simple Poisson distribution. A pre- 
liminary study of neural networks with a scale- free degree 
distribution showed that these networks demonstrate dy- 
namical properties qualitatively similar to properties of 
the networks studied above. 

Let us discuss possible experiments to test the pro- 
posed model. First, it would be interesting to ob- 
serve hysteresis and neural avalanches like those found 
in Sec. IIVI near the discontinuous phase transition. A 
similar discontinuous phase transition was revealed in ac- 
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tivations of living neural networks by a stimulus in recent 
works [m ■ Neural avalanches in these biological systems 
can be found by use of microelectrode arrays as described 
in [m, [3^. Second, the theory predicts that the emer- 
gence of global oscillations is a dynamical phase transi- 
tion. A strong enhancement of the response of a neural 
network to a periodic stimulus in the range of frequen- 
cies of these oscillations manifests this transition. These 
oscillations can be driven by an external stimulus, neu- 
ral pacemakers or noise. Though we demonstrated this 
behavior for ideal neurons, we believe that it is a uni- 
versal critical phenomenon if a sufficiently large group 
of neurons is involved in these oscillations. It would be 
interesting to observe experimentally this enhancement 
which in fact is stochastic resonance. This enhancement 
may be found, for example, in experiments similar to the 
experiments carried out by Fries et al. [48j who observed 
that neurons of macaque monkeys activated by the at- 
tendant stimulus show increased gamma-frequency (35 
- 90 Hz) oscillations. One can expect that a response 
of the neurons to periodic stimulus in the gamma-band 
frequency will be enhanced near critical attention above 
which gamma-frequency oscillations emerge. 

VIII. CONCLUSION 

In conclusion, based on experiments and ideas of cellu- 
lar automata we developed a model of noisy neural net- 
works with excitatory and inhibitory neurons and a com- 
plex network architecture. We considered neurons which 
are either inactive or fire a regular train of spikes with 
a given frequency (neurons with type 2 dynamical be- 
havior). In this model we took into account spontaneous 
neural activity, which plays the role of noise, the activa- 
tion of neurons by a stimulus, neural pacemakers, and 
interactions between neurons. We derived rate equations 
describing the evolution of the global neuronal activity. 
These equations are exact for infinite uncorrelated com- 
plex networks with arbitrary degree distributions, though 



for brevity, we presented results only for classical random 
graphs. This model has a complex phase diagram with 
self-organized active neural states, hybrid phase transi- 
tions, hysteresis phenomena and a rich array of behavior 
including decaying and stable oscillations, stochastic res- 
onance, and neural avalanches. We showed that global 
oscillations and stochastic resonance are intrinsic prop- 
erties of this non-linear dynamical system. The oscil- 
lations emerge when noise, i.e., the spontaneous neural 
activity, reaches a threshold level while stochastic reso- 
nance is a precursor of global oscillations. We also found 
that the network structure is important. The larger the 
connectivity the broader is the region with global oscil- 
lations. Our simulations revealed that even small groups 
of 50-1000 neurons display oscillations similar to large 
networks. 

Further development of the model can be done by tak- 
ing into account the real structure of neural networks 
(clustering, degree-degree correlations, modular struc- 
ture, and other structural properties), a dependence of 
firing rate on input, variability of synapses, evolution of 
network structure, for example, considering growing net- 
works, or variable strength of synapses, and so on. Apart 
from the perspectives discussed above, one can also apply 
this stochastic model to study communication between 
different groups of neurons and the transmission of ac- 
tivity from one group or module of neurons to another, 
taking into account noise and complex network architec- 
ture. 
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